rm(list=ls())
options(error=stop)

library(doParallel)
library(doRNG)
library(foreach)

seed = 1 # Random seed

source("2.0.0 MyFunc.R")

load("100716MSCM.RData")    # load dataset
data = data[,2:ncol(data)]  # delete intercept
pL0  = pL0 - 1
source("3.2 PMLE.R")
source("3.5 CallDR.R")

max.steps = c(5,10,20,50,100,200,500,1000)
para.name = c("theta", "phi", "gop", paste("eta",c("0*","1*")),
              "nll time")
results  = array(dim=c(1,length(max.steps), 6,pL0+3),
                 dimnames = list(1, max.steps, para.name,
                                 c("*=0", "*=1", "*=.", 1:pL0)))
names(dimnames(results)) = c("Simulation number","Max steps","Param",
                             "subscript")
for(i in 1:length(max.steps)){
        cat("step=",max.steps[i],"seed=",seed,"\n",sep=" ")
        mle = PMLE(data, count, K, seed, max.step = max.steps[i],
                    thres = 1e-4) 
        results[1,i,1, ]    = alpha.pmle = mle$theta.coef
        results[1,i,2, ]    = mle$phi.coef
        results[1,i,3, ]    = mle$gop.coef
        results[1,i,4:5, ]  = gamma.pmle =  mle$eta.coef
        results[1,i,6, 1:2] = c(mle$nll, mle$time)
        
        dr = DR
}    

save(results, file=paste("./RDataFiles/101016.maxsteps.seed",seed,".RData",
                         sep=""))

    
    